#include "../includes/need-a-name.hh"
#include "includes/initcollecter-constant-advection.hh"
#include <string>

InitCollecterConstantAdvection::InitCollecterConstantAdvection(){
  set_parameter("wid1",0.0);
  set_parameter("wid2",0.0);
  set_parameter("wid3",0.0);
  set_parameter("r",0.0);
  set_parameter("cen1",0.0);
  set_parameter("cen2",0.0);
  set_parameter("cen3",0.0);
  set_parameter("rectangle",0.0);
  set_parameter("circle",0.0);
}
InitCollecterConstantAdvection::~InitCollecterConstantAdvection(){
}

void InitCollecterConstantAdvection::initial_function(double* location, double pillarlen, int currentdim, int index){
  double componenttmp;
  double locationtmp[gridsetting::dimension];
  double basisInside[memory::dgorder*memory::intergnumber];      
  if(memory::dgflag==0){
    for(int tmpi=0; tmpi<memory::phynumber;tmpi++){
      if(para["rectangle"]==1.0){
	memory::variable[tmpi][index]=rectangle(location,tmpi);
      } else if(para["circle"]==1.0){
	memory::variable[tmpi][index]=circle(location,tmpi);
      }
    }
  }
  if(memory::dgflag==1){
    // std::cout<< "initial function for dg-algorithm" <<std::endl;
    basis.getBasis(basisInside);
    for(int tmpl=0; tmpl<gridsetting::dimension; tmpl++)
      locationtmp[tmpl]=location[tmpl];
    for(int tmpi=0; tmpi<memory::phynumber;tmpi++)
      for(int tmpj=0; tmpj<memory::dgorder;tmpj++){
    	componenttmp=0.0;
    	for(int tmpk=0; tmpk<memory::intergnumber; tmpk++){
    	  locationtmp[currentdim]=location[currentdim]+(-0.5+tmpk/(memory::intergnumber-1.0))*pillarlen;
    	  if(para["rectangle"]==1.0){
    	    componenttmp+=(2.0/memory::intergnumber)*(2.0*tmpj+1.0)/2.0*rectangle(locationtmp,tmpi)*basisInside[tmpj*memory::intergnumber+tmpk];	    
    	  } else if(para["circle"]==1.0){
    	    componenttmp+=(2.0/memory::intergnumber)*(2.0*tmpj+1.0)/2.0*circle(locationtmp,tmpi)*basisInside[tmpj*memory::intergnumber+tmpk];	    
    	  }
    	}
	memory::variable[tmpj*memory::phynumber+tmpi][index]=componenttmp;	
      }
    for(int tmpi=0; tmpi<memory::phynumber;tmpi++){
      if(para["rectangle"]==1.0){
    	memory::variable[memory::dgorder*memory::phynumber+tmpi][index]=rectangle(location,tmpi);
      } else if(para["circle"]==1.0){
    	memory::variable[memory::dgorder*memory::phynumber+tmpi][index]=circle(location,tmpi);
      }      
    }
  }
  // std::cout<< "moveing debug-point" <<std::endl;
}

double InitCollecterConstantAdvection::rectangle(double* location, int phyi){
  double cen[gridsetting::dimension];
  double wid[gridsetting::dimension];
  double tmp=0.0;
  if(gridsetting::dimension>=1){
    cen[0]=para["cen1"];
    wid[0]=para["wid1"];
  }
  if(gridsetting::dimension>=2){
    cen[1]=para["cen2"];
    wid[1]=para["wid2"];
  }
  if(gridsetting::dimension==3){
    cen[2]=para["cen3"];
    wid[2]=para["wid3"];
  }
  if(phyi==0){
    tmp=1.0;
    for(int tmpi=0; tmpi<gridsetting::dimension;tmpi++)
      tmp=tmp*(std::abs(location[tmpi]-cen[tmpi])<wid[tmpi]? 1.0 : 0.0);
  }
  if(phyi!=0)
    tmp=0.5;  
  return tmp;
};

double InitCollecterConstantAdvection::circle(double* location, int phyi){
  double cen[gridsetting::dimension];
  double radius;
  double tmp=0.0;
  if(gridsetting::dimension>=1)
    cen[0]=para["cen1"];
  if(gridsetting::dimension>=2)
    cen[1]=para["cen2"];
  if(gridsetting::dimension==3)
    cen[2]=para["cen3"];
  radius=para["r"];
  double distance=0.0;
  for(int tmpi=0; tmpi<gridsetting::dimension;tmpi++)
    distance=distance+std::pow(location[tmpi]-cen[tmpi],2);
  distance=std::sqrt(distance);
  if(phyi==0)
    tmp=(distance<radius? 1.0 : 0.0);
  if(phyi!=0)
    tmp=0.5;
  return tmp;
};
